function binv = custom_betainv(p, a, b)
% custom_betainv  - inverse function of the beta distribution
%
% FORMAT:       binv = custom_betainv(p, a, b)
%
% Input fields:
%
%       p           values for which the inverse beta CDF is computed
%       a, b        beta function parameters
%
% Output fields:
%
%       binv        inverse beta function for (p, a, b)

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 3 || ...
   ~isnumeric(p) || ...
   ~isnumeric(a) || ...
   ~isnumeric(b) || ...
    isempty(a) || ...
    isempty(b)
    error( ...
        'BVQXtools:BadArgument', ...
        'Missing or invalid argument given.' ...
    );
end
if isempty(p)
    binv = [];
    return;
end
osize = size(p);
if numel(p) == 1
    if numel(a) > 1
        osize = size(a);
    elseif numel(b) > 1
        osize = size(b);
    end
end
if ~isa(p, 'double')
    p = double(p(:));
else
    p = p(:);
end
if ~isa(a, 'double')
    a = double(a(:));
else
    a = a(:);
end
if ~isa(b, 'double')
    b = double(b(:));
else
    b = b(:);
end

% build output
binv = zeros(osize);

% special cases
binv(p < 0 | p > 1 | ~(a > 0) | ~(b > 0) | isnan(p)) = NaN;
binv(p == 1 & a > 0 & b > 0) = 1;

% find normal case
k = find(p > 0 & p < 1 & a > 0 & b > 0);
if ~isempty(k)
    if numel(a) > 1
        a = a(k);
    end
    if numel(b) > 1
        b = b(k);
    end
    if numel(a) ~= 1 || ...
        numel(b) ~= 1
    	y = a ./ (a + b);
    else
        y = a ./ (a + b) * ones(numel(k), 1);
    end
    p = p(k);
    
    % prepare loop
    y(y < eps) = sqrt(eps);
    y(y > (1 - eps)) = 1 - sqrt(eps);

    y_old = y;
    for lc = 1:10000
        h = (custom_betacdf(y_old, a, b) - p) ./ custom_betapdf(y_old, a, b);
        y_new = y_old - h;
        ind = find(y_new <= eps);
        if ~isempty(ind)
            y_new(ind) = y_old(ind) ./ 10;
        end
        ind = find(y_new >= (1 - eps));
        if ~isempty(ind)
            y_new(ind) = 1 - (1 - y_old(ind)) ./ 10;
        end
        h = y_old - y_new;
        h0 = (h == 0);
        if any(h0)
            binv(k(h0)) = y_new(h0);
            if numel(a) > 1
                a(h0) = [];
            end
            if numel(b) > 1
                b(h0) = [];
            end
            h(h0) = [];
            k(h0) = [];
            p(h0) = [];
            y_new(h0) = [];
        end
        if isempty(h) || ...
           (max(abs(h)) < sqrt(eps))
            break;
        end
        y_old = y_new;
    end
    if ~isempty(k)
        binv(k) = y_new;
    end
end

% reshape
binv = reshape(binv, osize);
